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ABSTRACT 

The aims of this paper are threefold: to increase the level of awareness within the 
shock capturing community to the fact that many Godunov-type methods contain 
subtle flaws that can cause spurious solutions to be computed; to identify one mech- 
anism that might thwart attempts to produce very high resolution simulations; and 
to proffer a simple strategy for overcoming the specific failings of individual Riemann 
solvers. 
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1 Introduction 

Over recent years, a plethora of shock-capturing schemes have been developed for 
the Euler equations of gas dynamics. During this period, it has emerged that one 
of the more successful strategies for designing a shock capturing scheme is to follow 
Godunov’s[19] lead and utilize a classic initial-value problem known as a Riemann 
problem [3]. Godunov assumed that a flow solution could be represented by a series of 
piecewise constant states. Thus, the numerical representation closely approximates 
the true solution near discontinuities, and regions of smooth flow are reasonably well 
approximated by a series of step functions. He evolved this discretized flow solution 
by considering the nonlinear interactions between its component states. Viewed in 
isolation, each pair of neighbouring states constitutes a Riemann problem, the solution 
to which may be found exactly [8]. The results from these separate Riemann problems 
may then be averaged so as to advance the flow solution through some time increment. 
Because it mimics much of the relevant physics, Godunov’s scheme results in an 
accurate and well-behaved treatment of shock waves. 

Although it provides the bedrock upon which most modern schemes are built, in 
its original form Godunov’s method is of limited use. Firstly, the scheme proves to 
be highly dissipative and so it requires an inordinately fine mesh to resolve complex 
shock-on-shock interactions. Secondly, since a Riemann problem has no closed form 
solution and can only be solved by some iterative method, Godunov’s scheme is 
significantly more expensive than schemes which employ ordinary finite-difference 
operators. 

One of the first people to address this second shortcoming was Roe[19], He argued 
that since the Riemann problems associated with Godunov’s method arise from an 
approximation of the data, it might be sufficient to find only approximate solutions 
to these Riemann problems, provided that they still describe important, nonlinear 
behaviour. His motivation being, that approximate solutions can be computed much 
more cheaply than exact solutions. Thus the industry of designing approximate Rie- 
mann solvers was born[14, 15, 4, 23]. Now whilst Godunov-type schemes are often 
held up to be models of robustness they can on occasions fail quite spectacularly. For 
example, when computing shock reflection problems, Roe’s method can sometimes go 
awry by admitting solutions for which the Mach stem is inexplicably kinked. The ex- 
istence of such failings partly explains why no consensus of opinion has been reached 
concerning the ideal Riemann solver. Whenever a new failing is unearthed, it adds 
fuel to the great Riemann solver debate; method X is better than method Y , because 
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of reasons A , B and C. It is our contention that, for all the current crop of Riemann 
solvers, at least one set of circumstances may be found for which the solver is found 
wanting; some failings are just more visible than others. 

In Section 2 we catalogue a number of situations in which anomalous behaviour is 
known to occur. This catalogue should serve to increase the general level of awareness 
within the shock capturing community to the current limitations of Riemann solver 
technology. At present, this awareness is not all that it should be, there are many 
instances in the literature where suspect numerical results are presented with either 
little or no adverse comment. We believe that one of the failings listed in our catalogue 
has hitherto gone unreported. In Section 3 we proffer a diagnosis of the mechanism 
which causes this new failing. 

At this juncture it should be noted that any foibles that a specific Riemann solver 
might have, may usually be controlled by the judicious use of a small amount of 
artificial dissipation. Indeed, it is worth pointing out that Woodward and Collela’s 
PPM scheme[2], which has proved itself to be more robust than most Godunov- type 
schemes, does in fact employ an elaborate artificial dissipation mechanism to supple- 
ment the dissipation provided via upwinding. As will be described in Section 4, we 
favour a strategy whereby the weaknesses of any one solver are overcome by combin- 
ing it with one or more complementary solvers. The main advantages of this approach 
over that of adding artificial dissipation are twofold. Firstly, it does not degrade the 
resolution of the base Riemann solver; it is possible to control certain instabilities by 
changing the flavour of the dissipation mechanism rather than increasing the absolute 
level of dissipation. Secondly, it does not necessitate a host of tunable parameters, 
and so this synergetic strategy does not negate the principal advantage of Godunov- 
type schemes over other shock capturing methods. Of course, we are left with the 
difficulty of deciding when to use one Riemann solver in preference to another; how- 
ever, we present a number of computations which suggest that this difficulty is not 
particularly bothersome. 

Finally, in Section 5, we list the main conclusions that we have drawn from this 
work. Note that in this paper we do not address the first shortcoming of Godunov’s 
method, namely its low resolution. Following van Leer[26], it is assumed that a high- 
order extension to a first-order method can always be achieved by pre-processing the 
data supplied to the Riemann solver. 
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2 A Catalogue of Failings 

We now present several instances where various Riemann solvers are known to give 
unreliable results. While most of these cases will be known to the aficionado, we 
believe that one of the cases which we are about to describe has hitherto gone un- 
reported in the literature. Although our catalogue is not exhaustive, we hope that 
it might save some investigators from the harrowing experience of spending weeks or 
even months searching for coding errors that simply do not exist. 

All the computational results shown in this section are for first-order schemes, 
since such methods have low resolution our calculations employed relatively fine 
meshes; for clarity, grids are drawn using every other fourth grid line. 

2.1 Expansion Shocks 

By far the most widely investigated failing is that some Riemann solvers do not 
satisfy an “entropy condition” such schemes can admit non-physical solutions such 
as expansion shocks. 0sher[14] has found a general condition for a scheme to be 
entropy satisfying when applied to scalar equations and he designates such schemes E- 
scheme 6 . At present, however, any extension to a system of equations contains a large 
amount of empiricism and must therefore remain suspect. Indeed, Godunov’s method 
is classified as an E-scheme but, as observed by Woodward and Collela[27], it can give 
rise to nearly discontinuous expansion fans near sonic points. The density contours 
shown in Figure 1 illustrate this deficiency of Godunov’s method quite clearly. These 
results are taken from the diffraction of a strong shock wave, M s = 5.09 with 7 = 1.4, 
around a 90° corner. 

In its basic form, Roe’s scheme is another solver that admits expansion shocks, 
however, several fixes have been proffered which cure Roe’s scheme of this particu- 
lar affliction[20, 5, 28]. Such fixes are typical of the way in which Riemann-solver 
deficiencies have been countered up to now. Whilst this strategy has proved reason- 
ably successful, it has a number of drawbacks. Sometimes a fix uses a parameter 
which must be retuned between problems and hence one of the major advantages of 
Riemann-based schemes over say artificial-dissipation schemes is lost. Alternatively, 
a scheme may require more than one fix, and it may be unclear how the different fixes 
interact with one another. 
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Figure 1: A strong shock diffracting around a corner gives rise to an expansion shock: 
(a) Density contours, (b) Computational grid. 

2.2 Negative Internal Energies 

Another situation, in which some Riemann solvers are found wanting, occurs when- 
ever the dominant energy mode is kinetic rather than thermal. For such solvers, 
the kinetic energy computed from a numerical approximation to the conservation 
laws of mass and momentum, can exceed the total energy computed via an approx- 
imation to the conservation law of energy. Thus they can yield negative internal 
energies, and hence negative pressures, which cause the scheme to fail. Einfeldt et 
al\ 5] call any scheme which can be guaranteed not to yield negative pressures, “pos- 
itively conservative”. They have shown that while Godunov’s scheme is “positively 
conservative”, the reverse is true for any Godunov-type scheme based on a linearized 
Riemann solver. Indeed, the basic form of Roe’s scheme is unable to cope with the 
test problem shown in Figure 1; the strength of the diffracting shock is sufficient to 
cause a negative pressure to be computed near the apex of the corner. Roe’s scheme 
may be made “positively conservative” by modifying its wave speeds, in essence, the 
scheme is made more dissipative by increasing the spread in velocity between the two 
acoustic waves [5]. 
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2.3 Slowly Moving Shocks 

Since shock capturing schemes do not resolve the internal structure of a shock wave, 
no physical significance can be attached to the discrete shock structure produced by a 
numerical scheme; methods are built upon the premise that shock profiles are mono- 
tone, the precise structure comes about as a matter of course and is not preordained. 
Unfortunately, Roberts has shown that the nature of the shock structure produced 
by a particular scheme can have a large bearing on how well the scheme copes with 
slowly moving shock waves[21]. Godunov-type methods fare quite badly in this re- 
spect, as the shock moves relative to the mesh, the shock profile flexes, perturbing 
the supposedly passive characteristic fields as it does so. 

Figure 2 shows a snapshot of the shock profile produced by Einfeldt's HLLE 
scheme[4], taken from the simulation of a shock wave which is moving slowly from 
left to right; the pre-shock state, (density, velocity, pressure), is (1,— 3.44,1), and 
the post-shock state is (3.86, -0.81, 10.33). Note that for a Courant number of one, 
it takes 50 time steps for this shock to traverse one mesh cell. The low frequency 
perturbations observed in this figure are produced to a greater or lesser extent by 
any scheme which attempts to “recognize” a shock wave. For fast moving shocks, 
the post-shock noise will be of a much shorter wavelength than is the case here, and 
will be effectively damped by the dissipation of the scheme. Roberts reports that 
Osher’s scheme[13] does not produce low frequency noise for slowly moving shocks, 
since it never connects two adjacent states by a shock, and he concludes that there 
may be advantages to using flux formulas that do not recognize the analytic shock 
jump conditions. 

Another situation, where the perturbation of a shock from its preferred profile 
results in perturbations on the passive characteristic fields, occurs whenever a shock 
crosses a discontinuity in mesh spacing[17]. But, in this case, sizeable perturbations 
may occur whatever the speed of the shock. 



Figure 2: Low frequency, post-shock oscillations occur for slowly moving shock waves. 
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2.4 The Carbuncle Phenomenon 

Several authors have by now reported a failing of Roe’s scheme which has been dubbed 
the “carbuncle phenomenon” [16, 12, 11]. For steady-state, blunt-body calculations, 
Roe’s scheme sometimes admits a spurious solution in which a protuberance grows 
ahead of the bow shock, along the stagnation line. It appears that this effect is 
more pronounced the more closely the grid is aligned to the bow shock. Also, a 
carbuncle is more likely to appear for high Mach number flows than for low Mach 
number flows. Figure 3 shows such a spurious solution, here the freestream Mach 
number was taken to be 10. Note that along the stagnation line, the bow shock is 
almost perfectly aligned with the grid. Consequently, parallel to the shock, Roe’s 
scheme will not add any dissipation via the contact and shear waves, to counteract 
perturbations that appear through the acoustic waves; this appears to be a recurring 
theme whenever Roe’s method fails. It is interesting to note that if Harten’s entropy 
fix[28] is applied to the contact and shear waves, any shortcoming of Roe’s scheme is 
invariably cured. However, there is no justification, either physical or mathematical, 
for applying this fix to these waves, it is just a convenient method for introducing an 
amount of artificial dissipation into the scheme. 

2.5 Kinked Mach Stems 

During the course of developing a mesh adaption scheme, we encountered a failing of 
Roe’s scheme which is not dissimilar to the “carbuncle phenomenon” [17]. When the 
reflection of a plane shock wave from a ramp lies in the double Mach reflection regime, 
the principal Mach stem is sometimes inexplicably kinked. Figure 4 shows a snapshot 
of the pressure contours taken during the reflection of a plane shock, M s = 5.5 with 
7 = 1.4, from a 30° ramp. The principal Mach stem is so severely kinked that it has 
given rise to a spurious triple point. Similarly strange results have been produced by 
Sawada[22], and by Itoh and Takayama[9]. As before, because of the way the Mach 
stem is aligned with the grid, there is probably insufficient dissipation added via the 
contact and shear waves to counteract perturbations that appear via the acoustic 


waves. 
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2.6 Odd-Even Decoupling 

By far the most insidious failing that we have come across has, we believe, gone 
unreported in the literature. During the course of producing very high resolution 
simulations, we have noticed a tendency for odd-even decoupling to occur along the 
length of planar shocks which are aligned with the mesh, for an example see Figure 1 1 . 
Of the Riemann solvers that we have at our disposal, this failing afflicts: an exact 
solver[24], Roe’s solver and Toro’s linearized solver[23]. We emphasize the fact that 
this phenomena only becomes apparent for very high resolution simulations which 
suffer some systematic pertubation. However, as will be shown below, the required 
perturbation can arise quite innocuously, so we suspect that this failing will prove 
fairly widespread once very high resolution simulations become commonplace due to 
increases in computer power. 

Now since we obtain our high grid resolution by means of a fairly complex mesh 
adaption scheme[17], it seemed reasonable to suppose that this odd-even decoupling 
was attributable to some coding error, but an exhaustive search for such an error 
proved fruitless. Subsequently, we have managed to reproduce this failing in a more 
controlled manner, as has a colleague using an independent code[10], so we have little 
doubt that this tendency for odd-even decoupling to occur, constitutes a genuine 
failing, rather than being the manifestation of some deficiency of our code. That 
said, our adaptive mesh scheme clearly exasperates this failing. In the next section 
we shall present a possible diagnosis of the mechanism which causes this mode of 
failure, here we merely present the evidence that it exists. 

Figure 5 shows several snapshots of the density contours from the simulation of 
a plane shock wave, M s = 6 with 7 = 1.4, propagating down a duct. For this 
calculation we have used Roe’s scheme, together with a nominally uniform grid of 20 
by 800 cells, with unit spacing, the centre-line of which is perturbed from that of a 
perfectly uniform mesh in the following manner 

y%jmid = Yjmid + 10“ 6 for i even, 

^ i,jmid — ^ jmid 10 for l odd. 

This perturbation to the grid centre-line promotes odd-even decoupling along the 
length of the shock. Note that the shock has propagated some 15 channel widths 
before the decoupling first becomes apparent, see frame (b). For this point in the 
calculation, Figure 6 shows a series of slices across the duct, for both the density and 
pressure fields, as one moves from the head to the foot of the shock. Interestingly, 
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within the shock the decoupling of the pressure field is out of phase with the decou- 
pling of the density field. As the shock continues to propagate down the duct, so 
the decoupling becomes progressively worse until the shock breaks down completely. 
However, at no stage in the calculation does the code blow up in the sense that it 
generates a floating point exception; it simply goes astray. 



(a) X s « 270 


(b) « 300 




(e) X s as 420 


(f) as 480 


Figure 5: Odd-even decoupling occurs for a shock propagating down a duct. 
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Figure 6: Transverse slices of the pressure and density fields for the shock shown in 
frame (b) of Figure 5. 
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3 Odd-Even Decoupling - A Diagnosis? 


As yet, the tools do not exist that would allow us to perform a rigorous nonlinear 
stability analysis for some shock capturing scheme applied to the Euler equations. 
However, it is possible to examine the way in which a scheme evolves certain sets 
of prescribed data, so as to ascertain its likely stability characteristics. Since some 
Riemann solvers allow odd-even decoupling to develop along the length of a plane 
shock, it might prove fruitful to examine how different schemes evolve sawtooth- type 
initial data. To this end, we consider the one-dimensional Euler equations with a 
passive component of shear velocity 



P ) 


' pv ' 

d 

pu 

d 

+ ■=- 

puv 

dt 

pv 

\ E ) 

By 

pv 2 + p 

K (E + p)v j 


The quantities p, p, u, v and E are density, pressure, the passive shear component 
of velocity, the velocity in the y direction, and the total energy per unit volume, 
respectively. For a perfect gas 


E = 


P 1 + yp { u2 + 

7-1 l 


where 7 is the ratio of specific heats. We assume that the computational mesh is 
uniform, with mesh spacing Ay, and that the discrete solution at time t n is given by 


71 „ , n _ 1 *ti 

Pj = P + p , Pj - p + P . 


Uj = u, 


V] = 0 , 


(2) 


if j is even, and 


Pj = P - P > Pj = P~ P 


u J - u, 


V] - 0, 


if j is odd. Here p n and p n are the amplitudes of the sawtooth profiles for the density 
and the pressure fields. We shall consider two schemes which may be expressed in 
the form 


W" +1 





( 3 ) 


where W is the conserved variable vector (p, pu, pn, E)\ and G n : is a first-order flux 

•J ' t 

function computed from the states W” and W” +1 . 
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3.1 Roe’s scheme 


Using Roe’s scheme[19], the interface flux for the system of equations (1) may be 
written 


G 


j+j 




where the wave speeds { A'*'}, wave strengths {cr*)} t and eigenvectors {e^} are given 

by 

A^ = v — a, A^ 2 ) = A* 3 ' = v, A^ 4 ) = u + a; 


(l) Ap — paAv 

a 1 ' = 

2a 2 


a (2) = A p-% 
a 1 


a 


( 3 ) 


pAu, 


1 Ap + paAv 

a K ' = 

2a 2 



1 1 


! ^ 


f°l 

1 

0 

w 


1 1 

e(U = 

u 

v — a 
yh — av y 

e (2) _ 
1 C 

It 

V 

^(u 2 + u 2 )J 

, e U) = 

, eU) — 

u 

v + a 
+ av j 


Here, quantities written as (•) are the so-called Roe averaged quantities, a and h are 
the Roe averaged sound speed and total enthalpy respectively, and A(t) represents 
the forward difference operator — (•) J . Now for our chosen data, 


(•)*! = Wy-l, and A(.) >+ . 


Therefore, 


k=4 

£ 

*= l 


«■% 


x(*) 

2 +J 


»<*> _ 


k=4 


2 *=i 


(*) 


x {k \ 


J-2 


Also, = Gj +1 , so the evolution scheme (3) may be written 


w™+i 

j 


At k = 4 
A v k = i 


-"7 + ££<*% 


v(fc) 
i+ 1 




1) 


which can simplified to 


w; +1 = w; + v y ^ 


( i \ 

U 

0 

\ ~ h ) 


( 4 ) 
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Where, u v is the Courant number . Recognizing that W" may be expressed as 
' y ' a 7 / 


0 \ 


' p n ) 

pu 

± 

up n 

0 


0 

£ J 


. p" 4. nifl" , 

h-i) + 2 P / 


and that by definition 

ip {t 2 

a 2 = ( 7 - i)(fc-y- T ), 

equation ( 4 ) may be manipulated to give, 


and 

P n+1 = p"(l — 2f/y). 

From which it can be seen that the initial perturbation to the pressure field is damped, 
provided that the CFL condition is met, that is 


Vy < 1. 

However, the form of the evolution for the density perturbation exposes a flaw r in 
Roe’s scheme; the density perturbation is fed directly from the pressure perturbation. 
Making the loose approximation that a remains constant, for a one-off disturbance 


p n = p° - ^p° [l + (1 - 2v y ) + (1 - 2 u y ) 2 + ■ • • + (1 - 2v y ) n ~\ 


Thus, 


r 



( 5 ) 


Therefore, if p° is of opposite sign to p° then for a one-off disturbance p will grow 
but it will remain bounded. But if the pressure field is continuously perturbed in 
a systematic manner, no matter how small the pressure perturbations, p will grow 
without bound, albeit slowly. For two-dimensional calculations, although we cannot 
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prove it, we suspect that a strong shock wave moving normal to the y direction 
provides this systematic perturbation. 

Firstly, it is interesting to note that the failing reported in Section 2.6 is only 
observed for strong shocks. For a strong shock wave, it seems reasonable to suppose 
that lfr| is more likely to be larger than |/5°|, than would be the case for weak shocks. 
So, even if p 0 and jP are initially of the same sign, they need not remain so, see (5). 
Now consider frame (b) of Figure 5 and the associated profiles shown in Figure 6, 
within the shock the odd-even decoupling of the pressure and density fields are indeed 
out of phase with one another, which is consistent with the observations made abo\e. 
Such behaviour will cause the local sound speed to vary along the length of the 
shock, its profile will exhibit a sawtooth perturbation which is in phase with that of 
the pressure field. Consequently, the individual segments of the shock will be moving 
alternately faster and slower than the nominal shock speed. Such movements will 
exaggerate the sawtooth perturbation to the pressure field along the length of the 
shock, but diminish that for the density field. The increased pressure perturbations 
will then promote an increase in the density perturbations as detailed above, and so 
the whole process repeats itself. 

Since there are two competing processes that affect the density perturbations, 
namely, the relative movements of the shock and the decoupling along the length of 
the shock, we cannot categorically state that Roe’s method is bound to break down. 
However, the weight of numerical evidence suggests, that at least for strong shocks 
Roe’s scheme will break down in the manner described here. Given our arguments, it 
should come as no surprise that Godunov’s method also exhibits a tendency to allow 
odd-even decoupling to occur along the length of a strong shock wave. Since it is 
the sweep parallel to the shock that primarily causes the instability, the differences 
between using an exact Riemann solver and Roe’s linearized solver for data that is 
nominally uniform should have little bearing on the growth of the instability. 

Finally, before moving on, it should be noted that none of the popular entropy 
fixes which are applied to Roe’s scheme cure this particular failing, excepting the 
case where Harten’s fix is applied to the shear and contact waves 2 ; simply altering 
the acoustic wave speeds can have no affect, because of the symmetry of the data, both 
waves will be changed by the same amount, and so the problem will persist. Also, 
moving to a high-order version of Roe’s scheme will not improve matters, because the 
odd-even decoupling will cause a high-order flux function to drop to the first-order 
function. 

2 To reiterate the comment made in Section 2.4, applying Harten’s entropy fix to the linearly 
degenerate wave fields has no mathematical or physical justification, it is merely a convenient way 
in which to add an amount of artificial dissipation. 
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3.2 Einfeldt’s HLLE scheme 

For Einfeldt’s HLLE scheme[4], the interface flux function is given by 


where 


& f+) G", , - &<"> G" 

pn Jt i 

J+s “ &(+) _ &(-) 


5rrv^( w " + . 


W”), 


6 (+) = max(0, A (4 ) ,,u ;+ i + a J+1 ), 6 ( > = max(0, a9\, v } - a,-), 

j > 2 ^ 2 

and A* 1 ), A^ 4 * are the two acoustic wave speeds from Roe’s method. Now for our 
chosen data, (») j+ l is equal to therefore 


&<t ) , = 


3+2 




and 


t ■ - -t\ 


-a, 


Using these signal weightings, it may be found that 

P n+1 =(l-2 


and 

P* +1 =(l-2;/ y )p", 


where, 


a&t 

Ay ' 


From which it can be seen that both the density and pressure perturbations are 
damped, provided that the CFL condition is satisfied. Just as important, however, is 
the fact that the pressure perturbation does not feed into the density perturbation, so 
we would not expect the HLLE solver to exhibit the odd-even decoupling that afflicts 
both Roe’s scheme and Godunov’s scheme; numerical experimentation confirms this 
expectation. 

It is our contention that any scheme for which it can be shown that the perturba- 
tion to the pressure field feeds the perturbation to the density field, will be afflicted 
by the odd-even decoupling shown in Figure 5. Thus it comes as no surprise to find 
that Toro’s linearized Riemann solver[23] is afflicted by this failing, but Liou’s ALISM 
scheme[ll] is not. The way is now open for some mathematician to perform a more 
rigorous analysis than we are able, so as to shed additional light on the mechanism 
which causes this particular failing. 
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4 An Adaptive Riemann Solver 


Having exposed many of the weaknesses of Riemann solvers, we now present a simple 
strategy, that we have found useful, for improving the robustness of Godunov-type 
schemes. In essence, we select the precise flavour of upwinding to match the local 
flow data, such that a particular Riemann solver is only employed in those situations 
where it is known to give reliable results. By recognizing the limitations of any one 
solver it is possible to reap its advantages without suffering its attendant failings. 

Our synergetic strategy has a number of attractions, not least of which is that some 
favoured solver need not be jettisoned simply because it, occasionally, fails. However, 
it does introduce the difficulty of how to decide when to use one Riemann solver in 
preference to another. But it has been our experience that this added difficulty is not 
particularly bothersome, for we tend to combine a single high resolution Riemann 
solver with just one or two other solvers that prove more reliable under conditions 
which are fairly well-defined, and so a set of ad hoc switching functions suffice. For 
example, some of the worst failings of Riemann solvers occur in the vicinity of strong 
shock waves. To overcome such failings we employ Einfeldt’s HLLE scheme. Now it 
makes little sense to chop and change the choice of Riemann solver used along the 
length of a shock wave, since to do so would inevitably perturb a planar shock front. 
Hence, we apply this particular Riemann solver throughout the immediate vicinity of 
a strong shock. Thus the HLLE switching function need only locate the position of a 
shock wave, but such functions already exist in the guise of mesh refinement, monitor 
functions, 

A simple test that identifies those cell interfaces which are in the vicinity of a 
strong shock is to check whether or not 


I Pr ~ Pl\ 
min(pi,p T ) 


( 6 ) 


where a is some threshold parameter which is problem dependent and p T and pi refer 
to the pressures which act on the interface. If this condition is met, the two cells 
separated by the interface are flagged as lying within a strong shock. So, when it 
comes to computing cell-interface fluxes, if the cells either side of an interface are 
both flagged as lying within a strong shock, the flux is computed using the HLLE 
solver. Note that since numerical shocks are invariably smeared over several mesh 
cells, it is worth locating shocks using a projection of the flow solution on a grid which 
is coarser than that used for the calculation. On such a grid a shock will appear much 
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less smeared, and so the left-hand side of the above switching function will be a fair 
indication of its strength. Once a set of cells have been flagged on this coarse mesh, 
the flags may be prolongated to the actual computational mesh so as to find those 
cells which lie in the vicinity of a shock wave. 

Before proceeding further, several observations should be made. Firstly, although 
the HLLE solver is adjudged to be a low resolution scheme, it does in fact resolve 
isolated shocks as well as an exact Riemann solver. Consequently, using this, robust 
solver in the vicinity of strong shock waves does not necessarily pollute a scheme’s 
resolution, as would be the case if artificial dissipation were used to augment the 
dissipation provided via upwinding. For many inviscid calculations the amount of 
pollution proves to be negligible, and whilst some degradation would be expected for 
the case of a strong shock interacting with a boundary layer, it may well be unneces- 
sary to employ the HLLE solver in such a situation because of the extra dissipation 
provided by the real viscous terms. Secondly, although the HLLE switching function 
requires a tunable parameter a, the retuning of this parameter is less involved than 
the retuning of an artificial dissipation mechanism; in general, it is far simpler to 
determine where extra dissipation should be added, than it is to determine how much 
extra dissipation to add. For many problems, assuming that shocks are located as 
described above, a sensible threshold on the shock strength can be specified a pri- 
ori. Lastly, it should be noted that our strategy of switching Riemann solvers may 
not prove suitable for those implicit schemes which require that the numerical flux 
function be differentiable. 

Figure 7 shows how the HLLE solver may be used to correct the tendency of 
Roe’s scheme to produce kinked Mach stems, c.f. Figure 4. For this calculation 
the HLLE switching function was tuned such that it would only be activated by the 
incident shock, and the principal Mach stem; o was simply set to half the strength of 
the incident shock. Note that apart from the region near the Mach stem, these new 
results are very similar to the old ones. This shows that the HLLE scheme has had no 
adverse affect on the the resolution of Roe’s scheme. Similarly, Figure 9 shows how the 
carbuncle phenomena may be circumvented, c.f. Figure 3. Here we have restricted the 
HLLE solver to cells near the stagnation line in order to demonstrate how localized the 
failing of Roe’s scheme really is. In practice, however, we would advocate using the 
HLLE scheme along the whole length of the bow shock, so as to maximize robustness 
without compromising resolution. Again, a sensible value of a can be found a priori 
by using some large fraction of the shock strength along the stagnation line which 
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can be estimated, given the freestream Mach number, by assuming the flow is locally 
one-dimensional. As shown in Figure 8, the HLLE solver may also be used to good 
effect to prevent Godunov’s scheme from admitting expansion shocks, c.f. Figure 1. 
Here we have employed the HLLE solver along the sonic line, and in regions where 
the expansion waves are strong. 

Having presented the gist of our strategy, we see little point in trying to sell 
a particular combination of solvers. Starting with some high resolution Riemann 
solver, whose choice will inevitably be a matter of personal taste, the correct com- 
bination of solvers will depend both on that schemes weaknesses and on the specific 
application in hand. In turn, the combination of Riemann solvers will dictate the 
choice of switching functions. Therefore, we shall resist the temptation to recom- 
mend a specific course of action, instead, we present two simulations that show how 
an adaptive Riemann solver might be used to good effect. Briefly, both simulations 
were done using the two-dimensional analogue of the one-dimensional Euler equations 
given in Section 3. These equations were integrated using the two step finite-volume 
scheme which is attributable to Hancock[25]. This scheme employs van Leer’s MLISCL 
approach [26] to achieve a second-order extension to Godunov’s method, hence differ- 
ent Riemann solvers may be slotted directly into the scheme so as to vary the flavour 
of the upwinding. Although the calculations were performed using an adaptive mesh 
algorithm[17, 18], the mesh refinement monitor function was such that the calcula- 
tions employed a nominally uniform cartesian mesh. 

Our first example concerns the simulation of a strong shock wave diffracting 
around a 90° corner, the shock Mach number and the ratio of specific heats are 
5.09 and 1.4 respectively. We have computed this test problem using a combination 
of three different Riemann solvers; Toro’s linearized Riemann solver was used to per- 
form the MUSCL reconstruction step of Hancock’s scheme as described by Quirk [18], 
and the upwinding step was performed by adaptively selecting between Roe’s solver 
and the HLLE solver. The parameter a used by the switching function (6) was set to 
1 so as to limit the HLLE solver to the incident and diffracted shock fronts, and to a 
small region near the apex of the corner. Figure 10 shows a Schlieren-type snapshot 
taken from this simulation, the different shades of grey depict the magnitude of the 
gradient of the density field, the darker the shade the larger the magnitude; details of 
this shading procedure are given in Appendix A. Here, it is not our intention to assess 
the accuracy of these results, the interested reader may do this using the experimental 
results of Bazhenova et al\ 1], and the computational results of Hillier [7] . Instead, we 
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wish to illustrate the fact that certain Riemann-solver failings, if left unaddressed, 
can place an upper limit on the resolution of simulations that may be performed. 

Consider the consequences of disabling the HLLE switching function, so that Roe’s 
solver alone is used for the upwinding stage of Hancock’s scheme. The tendenc} r of 
Roe’s solver to allow odd-even decoupling to occur along a planar shock wave which 
is aligned with the grid will, sooner or later, cause this simulation to come to grief, 
see Figure 11, thus precluding the possibility of performing very detailed simulations. 
By way of comparing the resolution of these two sets of results, for Figure 10 there 
are 560 mesh cells from the apex of the corner to the point where the Mach stem 
meets the wall, for Figure 11 there are only 120 cells. The question of whether or 
not our adaptive mesh algorithm contains some flaw which exasperates the odd-even 
decoupling is largely academic. The fact remains, that Roe’s solver is susceptible to 
this mode of failure whereas the HLLE solver is not. Whether the initial stimulus 
comes from a distorted mesh as in Section 2 or from some component of the mesh 
adaption scheme, as seems likely here, is immaterial. 

So as not to leave the impression that the above shortcoming is somehow pecu- 
liar to Roe’s method, we present a second set of results which are taken from the 
interaction of a planar shock wave with a half-diamond; the shock Mach number is 
2.85, the ratio of specific heats is 1.4, and the angle at the apex of the diamond is 
90°. As before, we have run this test problem using a combination of three different 
Riemann solvers, but this time we have substituted an exact Riemann solver[24] in 
place of Roe’s linearized solver. Figure 12 shows a Schlieren-type snapshot from this 
calculation, note that some 800 cells cover the width of the diamond so this calcula- 
tion is well resolved. Also, as an aside we note that the quality of these results may 
be gauged by comparing them with the experimental results given by Glass et ai[ 6]. 
Once again, if the HLLE switching function is disabled, the simulation is ruined by 
the odd-even decoupling that occurs along the length of the incident shock, see Fig- 
ure 13. Note that this second calculation is of lower resolution than the first, only 
400 cells cover the width of the diamond. 

In this section we have attempted to show that the robustness of Godunov-type 
schemes may be improved by employing an adaptive Riemann solver, where the spe- 
cific flavour of upwinding is altered to suit the local flow conditions. If used sensibly, 
this strategy can overcome most known failings of individual solvers. Despite our 
efforts, we recognize that the majority of shock-capturing practitioners will continue 
to use artificial dissipation as a band aid to fix a particular Riemann solver at the 
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first signs of any failing, simply because it is expedient to do so. Whilst we find this 
approach disappointing, the principal aim of this paper is to emphasize the fact that 
most of the Riemann solvers that are in common use must be augmented in some 
way, if they are to be used for the purpose of producing, genuinely, high resolution 
simulations of shock hydrodynamic phenomena. 


(a) 


(b) 


Figure 7: The HLLE scheme can be used to circumvent the tendency of Roe’s method 
to produce kinked Mach stems: (a) Pressure Contours, (b) HLLE switching function. 



(a) 


(b) 


Figure 8: The HLLE scheme can be used to prevent Godunov’s method from produc- 
ing expansion shocks (a) Density Contours, (b) HLLE switching function. 










Figure 10: A Schlieren-type snapshot from the diffraction of a strong shock wave 
around a 90° degree corner. 




Figure 11: On its own, Roe’s approximate Riemann solver cannot be used to repro- 
duce the resolution of the simulation shown in Figure 10. 
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Figure 12: A Schlieren-type snapshot from the interaction of a planar shock wave 
with a half-diamond. 



Figure 13: On its own, an exact Riemann solver cannot be used to reproduce the 
resolution of the simulation shown in Figure 12, 
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5 Conclusions 

Unless the dissipation provided via upwinding is augmented by some other mecha- 
nism, any Godunov-type scheme built upon a single Riemann solver will probably be 
flawed. For example, when subjected to some small but systematic form of pertur- 
bation, most popular Riemann solvers for the Euler equations, including the exact- 
solver, cannot prevent odd-even decoupling occurring along the length of a strong 
shock wave which is aligned with the computational mesh. Thus far, this flaw has 
gone largely unnoticed simply because it is only exposed by very high resolution sim- 
ulations. However, given that the required perturbations can arise quite innocuously, 
this mode of failure should prove fairly widespread once genuinely high resolution 
simulations become commonplace due to increases in computer power. 

Although most flaws can be controlled by the judicious use of a small amount 
of artificial dissipation, to do so necessarily leads to a reduction in the resolution of 
the scheme. We favour an alternative approach, whereby the the failings of any one 
Riemann solver are circumvented by combining it with one or more complementary 
solvers. In essence, we advocate selecting the precise flavour of upwinding to suit 
the flow data. Admittedly, this synergetic strategy is not as aesthetically pleasing as 
having a single Riemann solver for all occasions, but we have shown that it can be 
made to work quite effectively. Besides which, Riemann solvers are sometimes touted 
as being a solution-adaptive technique, so the concept of an adaptive Riemann solver 
is not that contrived. 

Looking to the future, it is to be hoped that genuinely multi-dimensional Rie- 
mann solvers will overcome many of the shortcomings of today’s dimensionally-split 
schemes. However, given the way in which the present shortcomings have been stum- 
bled across, these multi-dimensional schemes may themselves arrive complete with 
subtle failings with which to ensnare the unwary. 
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A Schlieren-type Plots 


The plots shown in Figures 10 and 12 depict the magnitude of the gradient of the 
density field, viz 


N 



+ 



2 


and hence they may be viewed as idealized Schlieren images. So as to accentuate 
weak flow features the following nonlinear shading function has been used, 

shade = exp(-kxp) 

Here, k is a constant and t/ms a weighting function given by 


ggj ~ jVgjo 
IW>|,-|V/>|o’ 


where 


|V/>|0 — ^0 

l^^ll = p\maxi 

ko and k\ being constants. Note that shade is limited to values between 0 and 1, 
so for a 24 bit colour graphics system the grey level shade may be converted to an 
< R,G,B > triplet using 

< 255 * shade , 255 * shade , 255 * shade > . 

For both Figures the constants k , & 0 , and ki were set to 5, 0.05 and —0.001 respec- 
tively. 


REPORT DOCUMENTATION PAGE 


jh. : rpoort.ng Durden +cr this :o'iection of information v to 

ihermg arc maintaining the Data needea. ana cc-motetma ano r e 
lection *» tnforrraticn/rri^aing sugge^tcr* Kt reducing this D-raen 
*15 Highway, Suite *204. ^rJmgtcr. ■/& 22202-4302, and to the 3t^*:e : 


1. AGENCY USE ONLY (Leave blank) 



Form Approved 
OMB Nc 0704-0188 


<eAirg tnitfjCTiO"s iearcntna e» st ng cata sources 
j.r.q th'«s bwraen estimate Of an, ::ner assert c’ 
r'o'Tiat on OwatiOh* -ana ^eoc'tv 'Z'b jetters^- 
ct (0?C4.:t8S) .Washington. r-C 2S5C3 


3. REPORT TYPE AND DATES COVERED 


4. TITLE AND SUBTITLE 

A CONTRIBUTION TO THE GREAT RIEMANN SOLVER DEBATE 


5. FUNDING NUMBERS 


C NAS1-18605 
C NAS1-19480 


6. AUTHOR(S) 

James J. Quirk 


WU 505-90-52-01 


7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 
Institute for Computer Applications in Science 
and Engineering 

Mail Stop 132C, NASA Langley Research Center 
Hampton, VA 23681-0001 


9. SPONSORING /MONITORING AGENCY NAME(S) AND ADDRE55{£5) 
National Aeronautics and Space Administration 
Langley Research Center 
Hampton, VA 23681-0001 


8. PERFORMING ORGANIZATION 
REPORT NUMBER 

ICASE Report No. 92-64 


10. SPONSORING /MONITORING 
AGENCY REPORT NUMBER 

NASA CR-191409 
ICASE Report No. 92-64 


11. SUPPLEMENTARY NOTES 

Langley Technical Monitor: Michael F. Card 

Final Report 


12a. DISTRIBUTION /AVAILABILITY STATEMENT 

Unclassified - Unlimited 


Submitted to International 
Journal for Numerical 
Methods in Fluids 


12b. DISTRIBUTION CODE 


Subject Category 02, 64 


13. ABSTRACT (Maximum 200 words) 

The aims of this paper are threefold: to increase the level of awareness within the 

shock capturing community to the fact that many Godunov-type methods contain subtle 
flaws that can cause spurious solutions to be computed; to identify one mechanism 
that might thwart attempts to produce very high reolution simulations; and to proffer 
a simple strategy for overcoming the specific failings of Individual Riemann solvers. 


14. SUBJECT TERM5 


15. NUMBER OF PAGES 

33 


Riemann solvers: shock waves; numerical artifices 16. price code 

A03 


17 SECURITY CLASSIFICATI ON 18. SECURITY CLASSIFICATION 19. SECURITY CLASSIFICATION 20. LIMITATION OF ABSTRACT 
* OF REPORT OF THIS PAGE OF ABSTRACT 

Unclassified Unclassified 


OF ABSTRACT 


NSN 7540-01-280-5500 


Standard Form 298 {Rev 

Prescribed by ANSI Sid Z39-18 
298-102 


2-89) 









